AD-A0B1  094 

UNCLASSIFIED 


OREGON  STATE  UNIV  CORVALLIS  SCHOOL  OF  OCEANOGRAPHY  F/ft  », 

BACKSCATTERINS  PROSRAMS  FOR  SPHERICAL  TARGETS. <U) 

APR  79  R  K  JOHNSON*  L  FLAX*  0  STANOLEY  N00014-76-C-0067 

REP-79-7  u. 


mm 


Reproduction  in  whole  or  in  part  is 
permitted  for  any  purpose  of 
the  United  States  Government 


Unclassified 


REPORT  DOCUMENTATION  PAGE 


SECURITY  CL  ASStFlC  A  TION  <jF  THIS  PAGE  '+h*n  Dm.*  EntmrtJ) 


HEAD  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 


r  .  .  (2.  GOVT  ACCESSION  NO.  3.  RECIPIENT  S  CATALOG  MjMitR 


s.  type  of  report  »  period  covered 


PERFORMING  ORG.  REPORT  number 


.ONTRACT  or  GRANT  NUMB£Rf*J 


J14-76-C-/)067| 


N  AM  £  AND  ADDRESS 


School  of  Oceanography  / 

Oregon  State  University 

Corvallis,  OR  97330  . 


II.  CONTROLLING  OFFICE  NAME  AND  ADDRESS  f 

Office  of  Naval  Research  (  /  ( 

Ocean  Science  £  Technology  Division  V 

Arlington,  VA  22217 


14.  MONITORING  AGENCY  NAME  &  ADDRESSfif  dHloront  from  Controlling  Otdcm)  I  15.  SECURITY  CLASS.  (ot  (him  import) 


16  DISTRIBUTION  STATEMENT 


/3L 


Approved  for  public  release;  distribution  unlimited 


17.  DISTRIBUTION  STATEMENT  (ot  thm  mbmtrmct  *nf*r*d  In  Block  20,  It  dl/tmrmnt  from  Rmport) 


18.  supplementary  NOTES 


Report  dated  April  1979 


13.  KEY  wOROS  (Confirm*  on  rmvmrmm  mldm  If  n*c*»m*ry  mid  identify  by  biock  mtmbmr) 


70  A9STRACT  (Confirm*  on  r*v*rmm  mid*  It  n*c*»m*fy  *"o  Idmntlfy  by  block  numb  *r) 

The  programs  presented  compute  the  acoustical  reflectivity  of  a  sphere  in  a 
fluid  medium  and  differ  in  the  allowed  physical  properties  of  the  spheres. 

The  principal  outputs  of  the  programs  are  plots  of  reflectivity,  R2,  as  a 

function  of  size-frequency,  ka.  The  reflectivity  is  defined  as  R  =  . 

^  /4  r 

where  is  the  backscattering  cross-section  of  the  sphere,  and  az/4  is  the 
backscattering  cross-section  of  a  completely  reflecting  sphere  of  radius  a. 

The  variable  k  is  the  wavenumber  in  the  medium.  These  plots  may  be  considered 


DD  t'™n  1473  EDITION  OF  I  NOV  •»  IS  0150L27E 

S/N  0103-01*-  8(01  i 


Unclassified 

SECURITY  CLASSIFICATION  OF  THlLf  AOE  f*N»«  Dj<*  *xim4) 

'  a  1 2-  a.  6  ? 


Unclassified _ _ 

HHiTY  CL  P  tw  a  r  ION  Oe  THl  i  PA  Dafa  F.i  fr+j) _ _  _ 

as  dimensionless  representations  of  target  strength  vs.  frequency.  The  con¬ 
versions  are  TS  =  10  log  (R2a2/4)  and  f  =  cka/2ita  with  (a)  in  meters#  (c)  in 
meters  per  second  and  (f)  in  kHz.  This  relation  for  frequency  is  such  that 
the  product  of  frequency  in  kHz  and  radius  in  mm  is  240  when  ka  =  1.  These 
programs  have  been  used  primarily  for  low  contrast  cases  at  relatively  low 
values  of  ka  (less  than  10) .  Some  adjustments  to  the  tolerance  parameters 
may  be  necessary  for  other  cases. 


Unclassified 


security  classification  op  this  pao*p»*w»  a*t*  «ii»i« 


BACKSCATTERING  PROGRAMS  FOR 
SPHERICAL  TARGETS* 


Richard  K.  Johnson1 
Lawrence  Flax2 
David  Standi ey1 


Reference  79-7 
April  1979 


G.  Ross  Heath 
Dean 


Oregon  State  University,  School  of  Oceanography 
2Naval  Research  Laboratory  |  ~ 


KTIS  Gfli-.il  ^ 

EDO  I u?  [ 

U:;  a:.  I 

Ju::'  i-  • .  '  • 


♦Supported  by  the  Office  of  Naval  Research 


Ey _ 


Introduction 


/  -  * 


These  programs  compute  the  acoustical  reflectivity  of  a  sphere  in  a 
fluid  medium.  The  programs  differ  in  the  allowed  physical  properties  of  the 
spheres.  The  prinicipal  outputs  of  the  nrograms  are  plots  of  reflectivity, 
R^,  as  a  function  of  size-frequency,  Vx.  The  reflectivity  is  defined  as 


tr 


where  is  the  backscattering  cross-section  of  the  sphere,  and  l2/4  is 
the  backscattering  cross-section  of  a  completely  reflecting  sphere  of  radius 
a.  The  variable  k  is  the  wavenumber  in  the  medium. 


These  plots  may  be  considered  as  dimensionless  representations  of 
target  strength  vs.  frequency.  The  conversions  are 

[(-  c> 

TS  =  10  log  (rC2/4)  and 


T"  “ 


TS  =  10  log  (RV/4) 
f  =  cka/2*a  * 


with  a  in  meters,  c  in  meters  per  second  and  f  in  kHz.  This  relation  for 
frequency  is  such  that  the  product  of  frequency  in  kHz  and  radius  in  mm  Is 
240  when  ka  =  1. 

These  programs  have  been  used  primarily  for  low  contrast  cases  (with  g  and 
h  near  one)  at  relatively  low  values  of  ka  (less  than  10).  Some  adjustments 
to  the  tolerance  parameters  may  be  necessary  for  other  cases. 


Program  Variables 


Input  Parameters 


DRATIO 

PRATIO 

SRATIO 

Z 


Symbol  Meaning 

6C  compresslonal  attenuation  in  sphere  (dB/wavelength) 

B$  shear  attenuation  in  sphere  (dB/wavelength) 

g  density  of  sphere/density  of  medium 

h  compresslonal  speed  In  sphere/speed  In  medium 

s  shear  speed  in  sphere/speed  in  medium 

ka  size- frequency  parameter 


Output  Parameters 


reflectivity  or  form  function 
reflectivity  squared 


Hii* 


1 


Fluid  Sphere 


The  program  SPHERF  is  a  simplified  version  of  SPHERE.  It  calculates 
r2  for  a  sphere  which  differs  from  the  fluid  medium  only  in  density  and 
compressional  sound  speed. 
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PROGRAM  SPHERF 

c««*«****«*«*s«*«x*x«****»**x*xx*»*«*»******x****»****x****«««*»**»**«*»**«*** 

C  SPHERF 

C  COMPUTES  BACK  SCATTERING  FORM  FUNCTION  FOR  A  FLUID  SPHERE 
C  WRITTEN  BY  LARRY  FLAX.NRL 
C  MODIFIED  BY  R  K  JOHNSON  AND  D  STANDLEY 
C 

c***************x***«************<************«*********************y********x 

DIMENSION  ZZ< 1000 ) » IE< 10 ) . GG  < 1000 ) 

COMPLEX  Cl »  TEMP .SUM. F 

DATA  N0/2HN  /. D/0. 1/ .EPS/O . 0005/ 

C - PLOT  COMMON.  ETC. 

COMMON  /PLT/P < 1 4 ) . IROT , I S I ZE » NXCH , NYCH . XFMT . YFMT . XLAB . YL AB 
DIMENSION  XFMT ( 2  >  *  YFMT  < 2)  »  XLAB( 10 ) » YLAB ( 10>  r PLID(3 ) 

DIMENSION  DATL8LC10) »PRL ( 3>  »DRL ( 3 > 

EQUIVALENCE  CDaTLBL»PLID<3> ) 

DATA  PLID/'SPHE' . 'RF  '.'  '/ 

DATA  P/7. .8. .2. .2. .0. .0. . 0 . . 0. • 0 . .0 . . 1 . . 10 . . 1 . » 10 . / 

DATA  XFMT/' <F4. '. '1)  ' . / . YFMT/ ' <F5 . ' . ' 0)  ' / .NXCH/4/ . NYCH/5/ 

DATA  PENUP/1/.  PENDHN/O/ 

DATA  PRL/'H  -  '.2*'  '/.  DRL/'G  -  '.2*'  '/ 

DATA  XLAB/'  KA  '.9*'  '/ 

DATA  Y LAB/'  R  < ' » ' DB >  ' .8* '  '/ 

INTEGER  PENUP.  PENDWN 

C _ 

C  GET  ACOUSTIC  PROPERTIES  AND  KA  RANGE 

c_ 

1  CALL  DATMSG 

CALL  DTMSG ( DATLBL )  IFOR  PLOT 

TYPE  100 

ACCEPT  101.  DRATIO.  PRATIO  IG.H 
20  TYPE  102 

ACCEPT  101.  ZFROM.ZTO.ZSTEP 
IEND»<ZTO-ZFROM)/ZSTEP+l .5 
IF < IEND  .LE.  1000) GO  TO  30 
TYPE  110 
GO  TO  20 
30  TYPE  103 

ACCEPT  104.  ITYP 
IF ( ITYP.EQ.NO ) GO  TO  50 
TYPE  105 
50  CONTINUE 

CI»<0..1.) 

G2LM-1000.  ISET  G2L  MINIMUM  ARBITRARILY  HIGH 

C»*««*«  START  LOOP 

DO  800  IZ-l.IEND 

Z»FLOAT  < IZ-1 ) XZSTEP+ZFROM 

ZL-Z/PRATIO 

TE1-0. 

TE2-0. 

TEMP-<0..0.) 

CALL  SBESJ(Z> 0 >BJ.D. IE ( 1 ) ) 

CALL  SBES J < ZL . 0 . B JL .D.IE(2>> 

CALL  SBESY (Z.0.BY.IEC3) ) 

DO  6  K>1. SO 
L*K-1 

X«<2*L+1 )*<-l )**L 

CALL  SBESJ<  Z.K  >BJ1 . 0. XE<4  > )  ISPHERICAL  BESSEL  FUNCTIONS 
CALL  SBESJ<ZL.K.BJL1.D.IE(5>> 

CALL  SBESY <Z.K.BY1.IE(6) ) 

DO  80  ICK-1.6 
IF(IE(ICK).EQ.0)G0  TO  80 
TYPE  1 06 . I CK . IE ( I CK ) 

80  IE(ICK)  -  0 

■  .  'j! 
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0036 
0037 
0058 
,  0059 
0060 
0061 
0062 
0063 
0064 
0065 
0066 
0067 
0068 
0069 


!  0072 
|  0073 
0074 
0075 
0076 
;  0077 
0078 
i  0079 
0080 


B JP*B JBL/Z-8 J 1  'FIRST  DERIVATIVES 

BJPL»BJL*L/ZL-BJL1 

BYP-BY#L/Z-BY1 

E2»ZL*B JPL /8  JL 

AN-E2/DRATI0 

R»BJ*AN-Z*BJP 

S«8Y*AN-Z*BYP 

U»R*R+S*S 

SUM* <  X/U ) *  <  CI#R*R-R*S ) 

TEMP-SUM+TEMP 
T*REAL<TEMP > 

T1*AIMAG(TEMP> 

GE1*ABS<  (T-TEl )/T) 

QE2=ABS ( <  T 1 -TE2 >/Tl) 

.INCLUDE  MORE  MODES  UNTIL  CHANGE  IS  LESS  THAN  EPS 
IF (QE1 • LE<EPS<AND< QE2 . LE .EPS  >  GO  re  13 
TE1-T 
TE2*Tl 

BJ  *  BJ1  ! BESSEL  OUTPUT  FOR  NEXT  ITERATION 

BJL  *  BJL1 
BY  *  8Y1 
CONTINUE 

F  *  2.  *  TEMP  /  Z  !F ORM  FUNCTION 

FI -REAL <F> 

F2-AIMAG<F> 

.MOBULUS*SGRT ( BACKSCATTERING/GEOMETRIC  CROSS-SECTION) 
G2»F1*F1+F2*F2 
G*S0RT<G2) 

ZZ ( IZ) “AL0G10 ( Z ) 

G2L*10 • *ALOG10< G2  > 

IF< ITYP.NE  .NO)TYPE  107 . Z fL f G f G2L 
GG<IZ)*G2L 

IF(G2L.LT.G2LM)G2LN=G2L  'SEARCH  FOR  MINIMUM  VALUE  OF  G2L 

CONTINUE 


PLOTTING  ROUTINE! 


TYPE  10S.G2LM 
ACCEPT  101.  YS 
IF(YS.GE.O)GO  TO  1 
ZLQGS=ALQG10( ZFRON ) 

ZL0GE=AL0G10(ZT0 ) 

MIN  *  INK <SIGN(ABS(ZL0GS)+0.96fZL0GS) > > 
IF < MIN  .GT.  0 ) MIN  *  MIN  -  1  ’ 

XMIN  *  10.  **  MIN 

MAX  *  INT<(SIGN<ABS<ZL0GE>+0.96fZLQGF>  )) 
IF  (MAX  .LT.  0  >MAX  =  MAX  +  1 

XMAX  =  10.  #*  MAX 

P<1)  *  7.  ! X  LENGTH 

P<2)  *  8.  ! Y  LENGTH 

P<3)  *  2. 

P(4>  *  2. 

P<5)  *  XMIN  \ 

P(6>  *  XMAX  ' 

P<7>  *  YS  ! YMIN 

P<8>  »  YS  +  80.  fYMAX 

P<?)  »  P(5)  ! XO 

P<10)  =  YS  ! YO 

CALL  AXIS<3)  'DRAW  AXES 

ENCODE  <  6 . 109 1  DRL ( 2 ) ) DRAT 1 0 
ENC0DE<6» 109*PRL<  2) iPRATIO 
CALL  PLOTXY <P(5) fP<8> .PENUP fO) 

CALL  PLOT(S.IXfIY) 

IX  *  IX  +  200 


!  YMIN 
!  YMAX 
(XO 
!Y0 

•DRAW  AXES 


! GO  TO  (XMIN. YMAX) 
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0121 

0122 

0123 

0124 

0125 

0126 

0127 

0128 

0130 

600 

0131 

0132 

0133 

0134 

0135 

0136 

100 

0137 

101 

0138 

102 

0139 

103 

0140 

104 

0141 

105 

0142 

106 

0143 

107 

0144 

108 

0145 

109 

0146 

110 

0147 

* 


CALL  STRING !IX,IY,28,PLID,0»2) 

CALL  STRING! IX. IY-60, 10. DRL,0,2> 

CALL  STRING !  IX ,  I Y-120 / 10 ,PRL  ,0,2) 

CALL  PL0TXY!P!5) ,P!7> ,PENUP.O> 

00  600  1-1,  IEND 
YY  =  GG( I ) 

XXX  -  ZZ<I> 

IP<I  .EQ.  t ) CALL  PLOTXY < XXX *  YY,  PENUP,  0 
CALL  PLOTXY! XXX,  YY,  PENDUN,  0> 

CALL  PLWAIT 
CALL  PLOT! 2, I 885) 

CALL  PLOT! 3) 

CALL  PLWAIT 
GO  TO  1 

FORMAT!'  ENTER  DRATIO,  PRATIO  ...  CG,H3  • 
FORMAT  (3F10.4) 

FORMAT!'  ENTER  ZFROM,  ZTO,  ZSTEP  !  ',*) 
FORMAT!'  TYPE  RESULTS?! Y/N>  ',*) 

FORMAT ! A2 ) 

FORMAT  !/'  KA  MODE  MODULUS 

FORMAT!'  REQ  PREC  NOT  ACHIEVED.  ROUTINE  •' 
FORMAT !F10.3,I8,F14.4,FX1.I> 

FORMAT! 'OPLOTJGIVE  START  OF  Y  SCALE  (MAX 
FORMAT !F6. 3) 

FORMAT!'  TOO  MANY  INCREMENTS  ...  TRY  AGAIN 
END 


PAGE  003 
!  LABEL  THE  PLOT 


!G0  TO  FIRST  POINT 


,*> 

20L0G  '/> 

12,'  ER=»',I2> 

VALUE*' ,F5.1, ' )  :  ',») 

//) 


R  (DB> 


Elastic  Sphere 


The  program  SPHERE  calculates  R2  for  a  elastic  sphere  in  a  fluid  medium. 
The  original  version  of  this  program  was  written  by  Larry  Flax. 
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PROGRAM  SPHERE 

C*«*************************************************************************** 
C  SPHERE 

C  COMPUTES  BACK  SCATTERING  FORM  FUNCTION  FOR  AN  ELASTIC  SPHERE 

C  WRITTEN  BY  CARRY  FLAX.NRL 

C  MODIFIED  BY  R  K  JOHNSONr  T  KEFFERi  AND  D  STANDLEY 
C 

C****************»Y************************* ******* ********* ****************** 

DIMENSION  ZZ<500>  »IE<10>  ,GG < 300) , A ( 3 , 3 > 

COMPLEX  Cl t TEMP>  SUM* F 
REAL  LMlf  LP1 

DATA  N0/2HN  /» D/0 . 1/ « EPS/O . 00 1/ 

C - PLOT  COMMON,  ETC. 

COMMON  /PLT/P( 14 ) . IROT , ISIZE , NXCH.NYCH. XFMT . YFMT , XLAB . YLAB 
DIMENSION  XFMT (2)  ,  YFMT <  2 ) , XLAB < 10) , YLAB< 1 0 >  , PH D < 3 ) 

DIMENSION  DATLBL1 10) , PRL < 3 ) , DRL < 3 > , SRL < 3 > 

EQUIVALENCE  < DATLBL ,PLID < 3 >  ) 

DATA  PLID/'SPHE' , 'RE  ','  '/ 

DATA  P/7. ,8. .2. .2. .0. ,0. ,0.,0.,0.,0.,1.,10.,1.,10./ 

DATA  XFMT/' (F4. '.' 1 )  ' . / » YFMT/ ' ( F5 . ' . '0 )  ' / , NXCH/4/ , NYCH/3/ 

DATA  PENUP/1/.  PENDWN/O/ 

DATA  PRL/ 'H  »  '.2*'  '/.  DRL/ '6  ’  '.2*'  '/ 

DATA  XLAB/'  KA  '.9*'  /.SRL/'S  =  '.2*'  '/ 

DATA  YLAB/'  R  < ' . ' DB >  ' . 8* '  '/ 

INTEGER  PENUP.  PENDWN 

C _ 

C  GET  ACOUSTIC  PROPERTIES  AND  KA  RANGE 

C 

1  CALL  DATMSG 

CALL  DTMSG< DATLBL;  'FOR  PLOT 

TYPE  100 

ACCEPT  101 ,  DRATIO ,PRATIO,  SRATIO  !G,H,S 

20  TYPE  102 

ACCEPT  101.  ZFROM  v  ZTO. ZSTEP 
IEND=<  ZTO—ZFROM ) /ZSTEP+1 • 3 
IF < IEND  .LE.  500) GO  TO  30 
TYPE  110 
GO  TO  20 
30  TYPE  103 

ACCEPT  104,  ITYP 

IF  < ITYP.EQ. NO) GO  TO  40 

TYPE  105 

*0  C*SRAT 10**2/ ( PR AT IQ# *2— 2 . * SR AT 10**2 )  ! POISSON ' S  RATIO 

30  CONTINUE 

CI=»<0. ,  1 .  ) 

G2LM=1000.  ‘SET  G2L  MINIMUM  ARBITRARILY  HIGH 

C******  START  LOOP 

DO  800  IZ=1 « IEND 

Z*FLQAT ( I Z-l ) *ZSTEP+ ZFROM 

ZL=Z/PRATIO 

ZS-Z/SRATIO 

TE1=0. 

TE2=0. 

TEMP=<0. ,0. ) 

CALL  SBES J< Z , 0 , BJ , D , IE ( 1  > ) 

CALL  SBES J (ZL . 0 . BJL , D , IE ( 2 ) ) 

CALL  SBESY <Z,0»BY»IE<3) ) 

CALL  SBES J ( ZS , 0 . BJS. D»IE<4>  > 

DO  6  K=1 , 50 
L=K-1 

X»<2*L+1)*< 1-2*M0D<L,2) ) 

CALL  SBESJ< Z, K  >  BJ  L , D , IE ( 5 ) )  'SPHERICAL  BESSEL  FUNCTIONS 
CALL  SBES J ( ZL ,K , BJL 1fD,IE(6)) 

CALL  SBES J ( ZS,K,BJS1,D,IE(7)> 
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0055  CALL  SBESY(Z»K»BYlfIE<8>) 

0056  BO  80  ICK=lf8 

0057  IF( IE( ICK) . EQ .0)GQ  TO  80 

0059  TYPE  106>  ICK r  IE  <  ICK ) 

0060  80  IE (ICK)  =  0 

0061  BJP=BJ*L/Z-BJ1  ! FIRST  DERIVATIVES 

0062  BJPL=BJL*L/ZL-BJL1 

0063  BJFS=BJS*L/ZS-BJS1 

0064  BYP“BY*L/Z-BY1 

0065  ZS2“ZS*ZS  (SECOND  DERIVATIVES 

0066  LM1 “FLOAT  <  L#(L-1 ) ) 

0067  LP1 “FLOAT  <  LX ( L  +  l ) ) 

0068  BJPPL=(LM1/(ZLXZL>-1. >XBJL+2 . XB JL1/ZL 

0069  BJPPS=(LM1/ZS2-1.)XBJS+2.XBJS1/ZS 

0070  A  < 1 » 1 ) “( BJL-2 . XCX8 JPPL > /( 1 . +2. XC  > 

0071  A< 1 ,3)=-2.XLPlX(ZSXBJPS-BJS)/ZS2 

0072  A(2» 1 )=ZLXBJPL 

0073  A  (  2  •  3  ) =LP 1 XB  JS 

0074  A(3.1)=2.X(ZLXBJPL-BJL) 

0075  A(3»3)=ZS2#BJPPS+FL0AT( (L+2)X(L-1 >  >*8JS 

0076  El=A(2»l)XA(3t3)-A(3.1)XA(2f3> 

0077  E“A(1, 1)XA(3»3)-A(1»3)XA(3»1 > 

0078  E2“E1/E 

0079  AN“E2/DRATI0 

0080  R=BJXAN-ZXBJP 

0081  S«BYXAN-ZX8YP 

0082  U-RXR+SXS 

0083  3LH»(X/U)*(CI*R*R-R*S) 

0084  TEMP»SUM+TEMP 

0085  T“REAL  <  TEMP ) 

0086  T1“AIMAG ( TEMP  > 

0087  QE1“ABS  <<T-TE1)/T) 

0088  0E2“ABS< <T1-TE2)/T1) 

C . INCLUDE  MORE  MODES  UNTIL  CHANGE  IS  LESS  THAN  EPS 

0089  IF <  QE1 . LE • EPS . AND  <QE2 . LE»£PS  >G0  TO  13 

0091  TE1=T 

0092  TE2“T1 

0093  BJ  =  BJ1  (BESSEL  OUTPUT  FOR  NEXT  ITERATION 

0094  BJL  “  BJL1 

0095  BJS  “  BJS1 

0096  BY  »  BY1 

0097  6  CONTINUE 

0098  13  F  =  2.  *  TEMP  /  Z  (FORM  FUNCTION 

0099  F1“REAL(F) 

0100  F2=AIMAG(F) 

C . MODULUS=SORT ( BACNSCATTERING/GEOHETRIC  CROSS-SECTION) 

0101  G2=F1*F1+F2*F2 

0102  G=SQRT(G2) 

0103  ZZ< IZ)=AL0G10( Z ) 

0104  G2L“10.*AL0G10<G2) 

0105  IF< ITYP . NE . NO ) TYPE  1 07 » Z »L r G » G2L 

0107  GG< IZ)=G2L 

0108  IF<G2L.LT.G2LM)G2LM=G2L  (SEARCH  FOR  MINIMUM  VALUE  OF  G2L 

0110  800  CONTINUE 

C _ 

C  PLOTTING  ROUTINE! 

C 

0111  TYPE  108»G2LM 

0112  ACCEPT  101 »  YS 

0113  IF(YS.GE.O)GO  TO  1 

0115  ZLOGS”ALOG10<  ZFROM ) 

0116  ZL0GE“AL0G10<  ZTO) 

0117  MIN  =  INT  < (SI GN( ABS  <ZL0GS)+0.96rZL0GS)  ) ) 

0118  IF (MIN  ,GT.  0 ) MIN  =  MIN  -  1 
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0120 

0121 

0122 

0124 

0125 

0126 

0127 

0128 

0129 

0130 

0131 

0132 

0133 

0134 

0135 

0136 

0137 

0138 

0139 

0140 

0141 

0142 

0143 

0144 

0145 

0146 

0147 

0148 

0149 

0150 

0152 

0153 

0154 

0155 

0156 

0157 

0158 

0159 

0160 

0161 

0162 

0163 

0164 

0165 

0166 

0167 

0168 

0169 

« 


XMIN  =  10.  «*  MIN 

MAX  =  INT  < (SIGN ( ABS ( ZLOGE >+0.96. ZLOGE )  )  > 

IF ( MAX  .LT.  0)MAX  =  MAX  +  1 
XMAX  *  10.  «*  MAX 
P<1)  =  7. 

P<2)  -  8. 

P<3>  -  2. 

P<4)  »  2. 

P<5)  *  XMIN 
P<6)  *  XMAX 
P < 7)  »  YS 
P<8>  =  YS  +  80. 

P<9)  »  P( 5) 

P<10>  *  YS 
CALL  AXIS ( 3 ) 

ENCODE (6. 109. DRL( 2  > ) DRAT  10 
ENC0DE<6r 109. PRL<2) ) PRAT 10 
ENC0DEC6. 109.SRL  (2) JSRATIO 
CALL  PLOT X Y <P(5)»P<B)» PENUP » 0 )  ! GO  TO  < XMIN » YMAX ) 

CALL  PLOT <5>IX>IY) 

IX  *  IX  +  200 

CALL  STRING( IX.IY.28r PLID. 0.2)  ! LABEL  THE  PLOT 

CALL  STRING( IX » I Y-60 . 10. DRL .Or 2> 

CALL  STRING <  IX. I Y— 120. 10 »PRL i 0 . 2) 

CALL  STRING( IX • I Y-l 80 » 10 . SRL . 0 > 2 ) 

CALL  PLOTXY(P ( 5) »P<  7 ) » PENUP .0) 

599  DO  600  1*1,  IEND 

YY  *  GG( I  ) 

XXX  *  ZZ(I) 

IF < I  .EO.  1 ) CALL  PLOTXY <  XXX  >  YY.  PENUP.  0)  !G0  TO  FIRST  POINT 

600  CALL  PLOTXY (XXX.  YY.  PENDUN.  0> 

CALL  PLUAIT 

CALL  PLOT <2. 1885) 

CALL  PLOT (3) 

CALL  PLUAIT 
GO  TO  1 

100  FORMAT < '  ENTER  DRATIO,  PRATIO.  SRATIO  ...  CG.H.S3  :  ',*> 

101  FORMAT  (3F10.4) 

102  FORMAT < '  ENTER  ZFROM.  ZTO.  ZSTEP  1  '.*> 

103  FORMAT ( '  TYPE  RESULTS?<Y/N>  '»*> 

104  FORMAT <A2> 

105  FORMAT  </'  KA  MODE  MODULUS  20L0G  '/) 

106  FORMAT < '  REQ  PREC  NOT  ACHIEVED.  ROUTINE  #'.I2»'  ER=',I2> 

107  FORMAT  <F10.3>I8.F14.4,F11,1)’ 

108  FORMAT ( 'OPLOT! GIVE  START  OF  Y  SCALE  (MAX  VALUE* ' »F5 . 1 » ' )  :  '.* 

109  FORMAT <F6. 3) 

110  FORMAT ('  TOO  MANY  INCREMENTS  ...  TRY  AGAIN'//) 

END 


!X  LENGTH 
!Y  LENGTH 


!  YMIN 
!  YMAX 
1X0 
!Y0 

! DRAW  AXES 


10 


R  (0B> 


R  (OB) 


Viscoelastic  Sphere 


The  program  ABSPHR  calculates  R2  for  a  viscoelastic  (absorbing)  sphere 
in  a  fluid  medium.  The  original  version  of  this  program  was  supplied  by 
Tokahi  Hasegawa. 
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OOOS 
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0007 
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0010 

0011 

0012 

0013 

0014 

0015 

0014 

0017 
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0019 

0020 

0021 

0022 


0023 
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0043 

0044 


PROSRAH  ABiPMP 

C999aa9898*a8988****9*«*9****:M******a************tt*ttttt********* 
C  A»8PHP 

c 

C  AC0U8TXC  BACK  SCATTERING  FROM  AN  A§80RiING  SPHERE: 

C  REFLECTION  FORH  FUNCTION  VS  KA 
C 

C  THIS  PROORAH  CALCULATES  THE  REFLECTION  FORH  FUNCTION  FOR  AN 
C  ABSORB I HO  ELASTIC  SPHERE  BASES  ON  TAKAHI  HASEQAMA •  YOSHIKO 
C  KITAOAWA  AND  YUHIKO  UATANABE .  SOUND  REFLECTION  FROH  AN 
C  ABSORBING  SPHERE.  J.ACOUST.  SOC .  AH. 42.  1298-1300.  1977. 

C 

C  MODIFIED  BY  R  K  JOHNSON  AND  D  STANDLEY 

088888888888888888888888889888888888888888888888888888888888888X888 

COMPLEX  X1.X2. JN.JB. A. B.G.D.C. E. H.F.HN.PH.CN 
COMPLEX  D1C. D2C.TC > PSC > DENSC . X3 
REAL  K1A.  K2A 

DIMENSION  XX ( 300  > . YP  <  300 ) • JN ( 30 ) . JB ( SO ) ,SJ(SO> 

DIMENSION  8N(S0) .ALP (30) . BET( 30 > .DJN (SO) . DJB(SO) 

DATA  YES/ ' YES  '/.  EPS/O . 001/, PX/34. 3737/ 

C - PLOT  COMMON,  ETC. 

COMMON  /PLT/P< 14 ) .1R0T,ISI2E.NXCH.NYCH,XFMT.YFHT.XLAB,YLAB 
DIMENSION  XFMT <  2 ) , YFHT <  2  > . XL AB  < 10 ) , YLAB ( 10  > , PL  ID ( 3  > 

DIMENSION  DATLBL(IO) ,AB1L(3) ,AB2L(3> , PRL ( 3 > , DRL ( 3 > , SRL ( 3 ) 
EQUIVALENCE  (DATLBL.PL 10(3) ) 

DATA  PLID/'ABBP' » 'HR  '/ 

DATA  P/7. . S., 2. ,2. ,0. .0. ,0. ,0. ,0. ,0. ,1.,10.,1.,10./ 

DATA  XFHT/'(F4.', '1)  ' » / , YFMT / ' ( PS . ' » ' 0 >  ' /.NXCH/4/ , NYCH/3/ 

DATA  PENUP/1/,  PENDWN/O/ 

DATA  AB1L ( 1 )/ ' AB1“' / , AB2L ( 1 ) /' AB2* '/.DRL ( 1 ) / ' 0  »  '/ 

DATA  PRL ( 1 ) / ' H  -  ' / , SRL ( 1 ) / ' 8  -  '/ 

DATA  XLAB/ '  KA  '.98'  '/ 

DATA  YLAB/'  R  <'.'DB)  '.88'  '/ 

INTEGER  PENUP,  PENDWN 

1  CALL  DATN8Q 

CALL  DTHSO(DATLBL)  IFOR  PLOT 

C _ 

C  OET  OPTION 

C  1)  CALCULATE  AND  SAVE  RESULTS  ON  FILE  FOR  LATER  PLOTTING 

C  2)  CALCULATE  AND  PLOT  RESULTS 

C  3)  READ  RESULTS  FILE  AND  PLOT 

C 

TYPE  104 

ACCEPT  103.  I OPT 

IF ( IOPT  . EQ .  3)00  TO  410 

2  TYPE  902 

ACCEPT  100. AB1.AB2 
TYPE  904 

ACCEPT  100 .DRAT 10, PRAT 10, SRAT 10  IO.H.S 
CC1  ■  1.  /  PRATIO 
CC2  -  1.  /  SRAT 10 

9  TYPE  907 

ACCEPT  100.  ZFROM,  ZTO,  ZSTEP 

IF(Z8TEP  .EQ.  0. 0 ) ZSTEP  •  0.23 

IMAX  •  IFIX(  (ZTO  -  ZFROM)  /  ZSTEP  +  1.0001) 

ZFUMAX  .LE.  300)00  TO  10 
TYPE  114 
00  TO  9 

10  YM  •  1000.  ILAROC  MINIMUM 

IHDR  -  0 

ET  ■  SECNDS ( 0 . )  (ELAPSED  TIME 

C 984944  START  LOOP 

DO  30  I-l.IMAX 

XX(I)  •  ZFROM  ♦  ZSTEP  8  FLOAT (I-l) 


•T1BLE. 

A 


0046 
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AAOI  002 


0047 

0040 

004? 

0090 

0091 

0092 

0093 

0094 

0099 

0097 

009? 

0041 

0043 

0049 

0044 

0049 

004? 

0071 
0072 
0073  21 

0074 

0079  211 

0074  3 

0079 
007? 

0091 
0092 
0093  22 

0094 

0099  221 

0094  4 

0097 
0099  29 

009? 

OOfO 
00?  1 
00?2 
0©?3 
00?4 
0099 
0094 
0097 
0099 
009? 

0100 

0101 

0102 

0103 

0104 

0109 

0104 

0107 

0109 

0109 

0110 

0111 

0112 

0113 

0114 

0119 

0114 

0117 

0119 

Oil? 


X  *  xxu> 

K1A  -  CC1  •  X 
K2A  -  CC2  »  X 

AL1  -  -Atl  *  K1A  /  AX 
AL2  -  -AB2  *  K2A  /  AX 

XI  •  CMALX (K1A.AL1 > 

X2  ■  CMALX (K2A.AL2) 

NM  >  1A1X<  X  ♦  9.0) 

I A  <  X  .LE.  3.41INM  -  7 
IF < X  .LI.  2.20)NM  -  4 
IA(X  .LI.  1 .93 )NH  «  9 
IA ( X  .LI.  1.12)NM  -  4 
IA  ( X  .LI.  0 . 9 )  NM  -  3 
NNM  -  NM+  3 

i a < aimao ( xi )  .ca.  o.)ao  to  21 

CALL  C999J<X1 .  JN.  NNM.  UR) 

IA ( UR  .IQ.  0)00  tO  3 
TYAI  905.  UR.  NNM.  Xi 
00  TO  3 

CALL  8JBC9(NNM<  REAL (XI).  DJN) 

DO  211  1C  -  1.  NNM 

JN(IC)  •  CMALX (DJN ( 1C ) .  0.) 

IA (A1MA0IX2)  .10.  0.)00  TO  22 
CALL  C8B9J(X2.  JB »  NNM.  UR) 

IA ( ICR  .IQ.  0)00  TO  4 
TYAI  909.  UR.  NNM.  X2 
00  TO  4 

CALL  SJBE9(NNM.  RCAL(X2>.  DJB) 
DO  221  1 C«  1.  NNM 
JB(IC)  •  CMALX < D JR <  1C ) .  0.) 
CALL  9JRI9 (NNM »  X.  9J ) 

CALL  SNRC9(NNM.  X.  9N) 

AJJ  «  0. 


IX  19  KA 
!  COMARI9910NAL 
I (HEAR 


I COMALIX  COM.  WAVE  NUM9IR 
I COMALIX  SHEAR  WAVE  NUMfCR 


(NUMBER  OA  M0DE9  TO  CALCULATE 
IAA9TER  IA  NOT  CMALX 


ANN  -  0. 

DO  20  N-l.NM 

T  *  AL0AT<N  -  1) 

N1  «  N 

M2  •  N  ♦  1 

TC  *  CMALX (T .0. ) 

A  *  TC  *  JN(N1)  -  XI  *  JN(N2) 
t  *  A  -  JN(N1) 

01C  «  CMALX ( 1 . .0. ) 

D2C  «  CMALX (2. .0. ) 

0  •  <X2««2  /  D2C  -  TC  $  (TC  -  D1C))  *  JN(N1)  -  D2C  *  XI  *  JN(N2) 

A  *  A  /  • 

0*0/9 

C  •  D2C  *  TC  •  (TC  f  D1C)  «  JB(Nl) 

I  •  (D2C  *  TC  *  TC  -  X2««2  -  D2C)  *  JB(Nl)  t  D2C  *  X2  •  JB(N2) 

H  "02C*  TC*  (TC4D1C)  «  ((D1C-TC)  *  JB(Nl)  +  X2  *  JB(N2>) 

B  -  C  /  I 
I  •  H  /  I 

ABC  •  CMALX< .9.0. ) 

DCM9C  »  CMALX(DRAT10.0. ) 

A  •  A9C*  (A  -  B)  *  X2**2  /((D  -  I)  *DCN9C) 

AJ  ■  T  *  fJ(Nl)  -  X  *  9J(N2 ) 

AN  »  T  «  9N(N1 )  -  X  *  SN(N2> 

HN  •  CMALX(9J(NI ) .  -9N(N1>) 

AM  •  CMALX(AJ.  -AN) 

CN  •  (CHALX(AJ.O, )-A*CMALX (9J(N1 > .0. > )/(A*MN-AH) 

ALA(Ni)  -  RCAL(CN) 

KT(N1)  «  AIMAO(CN) 

NS  •  N  -  1 

YN  •  (2.  *  FLOAT (NS)  A  1.)  *  (>1.)**N9 

AJJ  •  AJJ  ♦  YN  *  BET (N>  (ADD  UA  REAL  AND  1 MAO I NARY  MODES 
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0120 

0121 

0122 

0123 

0123 

0126 

0127 

0129 

0130 

0131 

0132 

0133 

0134 

0136 

0137 

0139 

0141 

0143 

0144 

0143 

0146 

0147 

0148 

0149 

0130 

0131 


0132 

0133 

0134 

0133 

0136 

0137 

0138 

0139 

0160 

0161 

0162 

0163 

0164 

0163 


0166 

0167 

0168 

0169 

0170 

0171 

0173 

0174 

0173 

0177 

0178 

0179 

0180 

0181 

0182 

0183 


PNN  -  PNN  +  YN  *  ALP  !  N  > 

C . TEST  FOR  CONVERGENCE 

T1  -  ABS( YN  *  BET (N> )  /  PJJ 
T2  -  ABS<YN  *  ALP  !  N  >  >  /  PNN 

IF<ABS(T1)  .LT.  EPS  .AND.  ABS<T2>  .LT.  EPS)GO  TO  201 
20  CONTINUE 
N  -  N  -  1 

IF!NM+2  .GT.  NNN)GO  TO  200 
NM  -  NN  ♦  1 


200 

201 


C  ,  •  •  • 


291 


GO  TO  25 

TYPE  113,X,NNM,NM 

YP< I )  -  2,0  *  SQRT(PJJ**2  +  PNNM2)  /  X  ! CONVERT  TO  FORM  FUNCTION 
YPL  -  20,  *  AL0G10< YP< 1 ) ) 

IF(YPL  ,LT ,  YM> YH  ■  YPL 
.TYPEOUT  OPTION 
ISM  «  IPEEKC 177570) 

IF < ISW  .LT.  0)G0  TO  30 
IF ( ISW  .EQ.  0 ) GO  TO  291 
IF(IHDR  , EO,  0>TYPE  110 
IHOR  -  1 
DT  -  SECNOS !  £T  > 

TYPE  lllf  Xf  YP! I >  »  20.*ALOG10!YP!I)>,  NNM»NM»  XI, 

ET  »  SECNDSCO. ) 

GO  TO  30 
IMAX  =*  I 
GO  TO  31 
CONTINUE 
TYPE  102, 


30 
31 

C - 

C  CHECK  OPTIONS 
C 


! READ  CONSOLE  SWITCHES 
! NORMAL 

IEXIT  LOOP  IF  SWITCHES-0 
! TYPEOUT  IF  SWITCHES  + 

X2»  DT 


<XX<I>»  YP(I),20.*ALOG10(YP< I )>,!*!, IMAX) 


GO  T0< 400,430,410) , IOPT 
400  TYPE  106  ! OUTPUT  FILE 

CALL  ASSIGN ! 1 »  » - 1 » ' NEW  '  ) 

WRITE < 1 , 107) AB1 , AB2, ORATIO  »PRATIO»SRATIO, YM , IMAX 
WRITE< 1 , 108) !XX!  I )  »  YP!  I  )  » 1*1 » IMAX) 

CALL  CLOSE < 1 ) 

CALL  EXIT 

410  TYPE  112  ! INPUT  FILE 


CALL  ASSIGN! 1, »-l , 'RDO ' > 

READ! 1, 107 )AB1,AB2,DRATI0,PRATI0,SRATI0,YM, IMAX 
READ! 1 ,108) !XX< I ) ,YP! I ) ,1-1 , IMAX) 


CALL  CLOSE! 1) 
ZFROM  3  XX! 1) 
ZTO  *  XX! IMAX) 


C - 

C  PLOTTING  ROUTINE) 

C 

430  TYPE  109, YM 

ACCEPT  904,  YS 
ZLOGS-ALOG 1 0  <  ZFROM ) 

ZLOGE *AL0G10 ! ZTO ) 

MIN  -  I NT !! SIGN (ABS!ZLOGS )+0.96,ZL0GS ) ) > 
IF!MIN  ,GT .  0)MIN  =  MIN  -  1 
XMIN  *  10.  *»  MIN 

MAX  -  INT!!SI6N!ABS! ZLOGE )+0.96, ZLOGE ) ) > 
IF! MAX  .LT.  0 )MAX  -  MAX  +  1 


XMAX 

a* 

to.  **  MAX 

P !  1 ) 

• 

7. 

!  X 

LENGTH 

P!2) 

m 

8. 

fY 

LENGTH 

P!3) 

m 

2. 

P!4> 

m 

2. 

P!3> 

m 

XMIN 

P!6) 

m 

XMAX 

f 


i 
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0184 

0185 

0186 

0187 

0188 

0189 

0190 

0191 

0192 

0193 

0194 

0195 

0196 

0197 

0198 

0199 

0200 

0201 

0202 

0203 

0204 

0205 

0206 

0207 

0209 

600 

0210 

0211 

0212 

0213 

0214 

999 

0215 

100 

0216 

102 

0217 

103 

0218 

104 

0219 

105 

0220 

106 

0221 

107 

0222 

108 

0223 

109 

0224 

110 

0225 

111 

0226 

112 

0227 

113 

0228 

114 

0229 

902 

0230 

904 

0231 

905 

0232 

906 

0233 

907 

0234 

* 


P< 7)  -  YS  !  YMIN 

P<8)  -  YS  +  80.  ! YMAX 

P<9)  ■  P<3>  1X0 

P<10>  -  YS  IYO 

CALL  AXIS (3)  (DRAW  AXES 

ENCODED*  103 *AB1L: 2)  )AB1 

ENCODE  : 6  *  1 03  *  AB2L : 2 ) )AB2 

ENCODE  <  6 , 103  *  DRL ( 2 ) ) DRATIO 

ENCODE < 6. 1 03 .PRL < 2 > >PRATIO 

ENCODE <6*103*  SRL ( 2 ) JSRATIO 

CALL  PL0TXY:P:5) *P:8) *PENUP*0> 

CALL  PLOTCSt IX* IY> 

IX  *  IX  +  200 

CALL  STRING ( IX*IY*28*PLID»0*2> 

CALL  STRING ( IX * IY-60* 10* AB1L*  0*2) 


CALL  STRING: IX* IY-120 • 10* AB2L*0>2 > 
CALL  STRING ( IX* IY-180»10*DRL*0»2) 
CALL  STRING< IX  *  I Y-240  » 10* PRL  *0  *  2 ) 


!  GO  TO  ( XMIN  r  YHAX  > 

!  LABEL  THE  PLOT 


CALL  STRING ( IX* I Y-300* 10 *SRL* 0 *2) 

CALL  PL0TXY<P<5) *P(7) .PENUP.O) 

DO  600  1-1*  IHAX  ! PLOT  THE  POINTS 

YY  -  20.  *  ALOG10( YP< 1 ) > 

XXX  »  AL0G10<XX(I) ) 

IF ( I  ,EQ.  1 ) CALL  PLOTXY ( XXX .  YY*  PENUP*  0)  (GO  TO  FIRST  POINT 
CALL  PLOTXY < XXX*  YY*  PENDUN*  0> 

CALL  PLWAIT 
CALL  PLOT <2*1885 ) 

CALL  PLOT < 3) 

CALL  PLWAIT 
GO  TO  1 

FORMAT  <  8F 10.5) 

FORMAT ( /3 ( 5X * ' KA ' *  6X  * ' YP ' *6X*'20LOG' > »/ . 3 <2X *F7 . 3*F8 . 4 . 2X *F7 . 1 )  ) 
FORMAT <F6. 3) 

FORMAT < '  OPTION:  DCALCSSAVE  RESULTS*  2>CALC»PL0T*  3) READ  '* 

1  ' RESULTS tRLOT  :  '*•) 

format: 12) 

FORMAT: 'OOUTPUT  FILE  '*•) 

F0RMAT:6E1S.5*I6) 

F0RMAT:6E13.5) 

format: -oplot:give  start  of  y  scale  :max  value-' *F6.i* ' ) 

FORMAT :/3X. 'KA' >6X. ' YP' *4X» '20L0G' *3X» 'MODE' *SX* ' X1R ' * 7X* ' XI I ' * 

1  7X. 'X2R' *7X*'X2I'*3X* 'SECONDS'/) 

FORMAT :F7  <3*F8«4*F7.1*I4*13*2X*4:iPE10.3)*  2X  *0PF4 . 0 ) 

format:'  input  file  '*•> 

FORMAT:'  YP  DID  NOT  CONVERGE  AT  X  -  '.F5.2* 

1  '  NNM  =*  '*13*'  NM  <■  '*I3> 

format:'  too  many  increments  ...  try  again') 
format:'  enter  abi*  ab2  :  '*•> 

F0RMAT:2F10.5) 

format: 'OCSBSJ  ERROR  ♦'*I4*2X. 'M0DE-'*I3*2X*2E16.7> 

FORMAT:'  enter  DRATIO*  pratio*  SRATIO  ...  CG.H.S3  :  '*») 

FORMAT:'  ENTER  ZFRQHr  ZTO*  ZSTEP  S  '*•) 

END 
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Fig.  4  Output  from  ABSPHR  for  a  case  with  shear  attenuation 
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Real  Spherical  Bessel  Functions 

SBESJ 

This  procedure  evaluates  the  spherical  Bessel  function  of  the  first 
kind,  jn(x),  for  real  arguments.  The  method  used  for  computation  depends 
on  the  argument  and  the  order. 

For  orders  0  and  1,  the  exact  relations  are 
Jn(x)  =  sin  (x)/x 

0  (1) 
j-j(x)  =  sin  (x)/x2  -  cos(x)/x. 

For  higher  orders  and  small  arguments,  the  procedure  uses  the  series 
approximations 


This  approximation  is  acceptable  only  when 
x2/D  <  2n, 

where  D  is  the  required  accuracy  of  the  procedure. 

For  all  other  cases,  the  procedure  uses  the  recurrence  relation 

i„  ♦  i(x)  ’  l£Lri  VK)  -  Jn  -  i<*>-  (3) 

This  relation  can  be  used  either  to  ascend  to  higher  orders  or  to  descend 
to  lower  orders. 

Ascending  recurrence  Is  used  for  x  >  2n  since,  in  this  case,  the  error 
in  Jn (x)  will  not  be  Increased  in  Jn  +  ^x). 

For  the  remaining  case  (  x  <  2n),  the  procedure  uses  Miller's  device, 
which  is  a  decreasing  recurrence  technique.  This  method  uses  the  approximation 

A  A 

jm(x)  *  0  and  _  -j (x)  =  1  for  some  m  >  n,  then  uses  decreasing  recurrence 

A 

to  calculate  j^(x)  for  0  £  1  <  m  -  1.  A  scale  factor  P  Is  calculated  from 
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j0(x)/50(x) ,  where  jQ(x)  is  determined  from  equation  (1).  The  correct  value 
jn(x)  is  the  P*j  (x).  This  sequence  is  repeated  for  higher  values  of  m  until 
successive  values  for  jn(x)  differ  by  less  than  D. 

SBESY 

This  procedure  evaluates  the  spherical  Bessel  function  of  the  second 
kind  (also  called  the  spherical  Neumann  function),  yn(x)  or  nn(x),  for  real 
arguments.  The  method  is  ascending  recurrence  and  presents  no  computational 
problems. 

These  programs  were  written  by  Richard  Johnson  and  Thomas  Keffer. 
Bibliography 

Abramowitz,  M.  and  I. A.  Stegun,  1970.  Handbook  of  Mathematical  Functions 
(U.S.  Government  Printing  Office)  435-478. 
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SUBROUTINE  SBESJ 
PURPOSE 

COMPUTE  THE  J  SPHERICAL  BESSEL  FUNCTION  FOR  A  GIVEN  ARGUMENT 
AND  ORDER 

USAGE 

CALL  SBE SJ<X.N.BJ.D.l£R> 

DESCRIPTION  OF  PARAMETERS 

X  -THE  ARGUMENT  OF  THE  J  BESSEL  FUNCTION  DESIRED 
N  -THE  ORDER  OF  THE  J  BESSEL  FUNCTION  DESIRED 
BJ  -THE  RESULTANT  J  BESSEL  FUNCTION 
D  -REQUIRED  ACCURACY 
IER-RESULTANT  ERROR  CODE  UHERE 
IER*0  NO  ERROR 
IER-1  N  IS  NEGATIVE 
IER>2  X  IS  NEGATIVE  OR  ZERO 
IER-3  REQUIRED  ACCURACY  NOT  OBTAINED 

IER“4  RANGE  OF  N  COMPARED  TO  X  NOT  CORRECT  (SEE  REMARKS) 
REHARKS 

N  MUST  BE  GREATER  THAN  OR  EQUAL  TO  ZERO'  BUT  IT  MUST  BE 
LESS  THAN 

20+10BX-X**  2/3  FOR  X  LESS  THAN  OR  EQUAL  TO  15 
90+X/2  FOR  X  GREATER  THAN  15 

SUBROUTINES  AND  FUNCTION  SUBPROGRAMS  REQUIRED 
NONE 

METHOD 

RECURRENCE  RELATION  TECHNIQUE  DESCRIBED  BY  H.  A.  ANTOSIEUICZ 
IN  HANDBOOK  OF  MATHEMATICAL  FUNCTIONS »PP . 438-439 
AND  452-453. THIS  CODE  PATTERNED  AFTER  THAT  FOR 
BESJ. 


SUBROUTINE  SBESJtX.N. BJ»D» IER) 
BJ-.O 

IF  <N) 10.20.20 
10  IER»1 
RETURN 

20  IF <X)30>30f 31 

30  IER-2 
RETURN 

31  IF(X-15. >32.32. 34 

32  NTEST»20.  +  10.*X-X«  2./3. 

GO  TO  36 

34  NTEST-90.+X/2. 

36  IF (N-NTEST >40.38.38 
38  IER-4 
RETURN 
40  IER«0 
Nl-N+1 
BPREV-.O 

IF  N  IS  0  OR  1.  COMPUTE  DIRECTLY. 
IF  (N-l)  42.43.49 
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0020 

42 

BJ-SIN<X)/X 

0021 

RETURN 

0022 

43 

8J»SIN(X)/<X*X)-C0S<X)/X 

0023 

r» 

RETURN 

U 

c 

IF  X  IS  VERY  SHALL  USE  ASCENDING  SERIES 

0024 

A 

47 

IF ( XtX/D-FLOAT ( 2tN ) ) 300 1 ZO » 30 

L 

c 

IF  X  IS  LARGE »  USE  INCREASING  RECURRENCE 

0023 

L 

r 

SO 

IF<X-FL0AT<2*N>)60»210,210 

u 

c 

IF  X  IS  HIDDLINf  USE  DECREEING  RECURRENCE 

0026 

U 

60 

HA-X+8. 

0027 

NB*N+IFIX<X>/4+2 

0020 

MZERO-HAXO(HAfMB) 

L 

c 

P 

SET  UPPER  LIHIT  OF  H 

0027 

V. 

HHAX*NTEST 

0030 

DO  170  M»HZER0,MHAX>3 

0031 

FN1-1. 

0032 

FH«.0 

0033 

DO  160  K-l»M-i 

0034 

MK-H-K 

0033 

BHK»FL0AT<2*HK+1 )*FM1/X-FH 

0036 

FN-FM1 

0037 

FH1-BHK 

0038 

IF (NK-N— 1 ) 160* 140, 160 

0037 

140 

BJoBHK 

0040 

r* 

160 

CONTINUE 

t 

c 

p 

SCALE  FACTOR 

0041 

u 

P>SIN<X)/'(X*BHK) 

0042 

BJ»P*BJ 

0043 

IF  <  ABS  <  B J-BPREV ) -ABS  <  D*B J ) >200, 200 . 170 

0044 

170 

BPREV-BJ 

0043 

IER-3 

0046 

f* 

200 

RETURN 

u 

c 

p 

INCREASING  RECURRENCE 

0047 

L 

210 

BJA«SIN(X)/X 

0048 

B JB-SIN < X ) / <  X*X ) -COS  <  X  >  /X 

0047 

K»1 

0030 

220 

T«FLQAT(2*K+1)/X 

0031 

BJC-T4BJB-BJA 

0032 

K-K+l 

0033 

IF  (K-N)  230,240*230 

0034 

230 

BJA-BJB 

0035 

BJB-BJC 

0036 

GO  TO  220 

0037 

240 

BJ-BJC 

0038 

P 

RETURN 

c 

ASCENDING  SERIES: 

0037 

300 

BJ«1.0 

0060 

DO  380  IFAC-1,N 

0061 

BJ-BJ*X/( FLOAT (27IFAC+1 ) ) 

0062 

350 

CONTINUE 

0063 

RETURN 

0064 

END 
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0003 
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0005 
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0008 

0009 
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0011 

0012 


C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 
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c 

c 

c 


c 

c 

c 

c 

c 

c 


c 

c 
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SUBROUTINE  SBESY 
PURPOSE 

COMPUTE  THE  Y  SPHERICAL  BESSEL  FUNCTION  FOR  A  GIVEN 
ARGUMENT  AND  ORDER 

USAGE 

CALL  BESY ( X .N. BY . IER  > 

DESCRIPTION  OF  PARAMETERS 

X  -THE  ARGUMENT  OF  THE  Y  BESSEL  FUNCTION  DESIRED 
N  -THE  ORDER  OF  THE  Y  BESSEL  FUNCTION  DESIRED 
BY  -THE  RESULTANT  Y  BESSEL  FUNCTION 
IER-RESULTANT  ERROR  CODE  WHERE 
IER=0  NO  ERROR 
IER=1  N  IS  NEGATIVE 
IER-2  X  IS  NEGATIVE  OR  ZERO 
IER-3  BY  HAS  EXCEEDED  MAGNITUDE  OF  10**36 

REMARKS 

X  MUST  BE  GREATER  THAN  ZERO 
N  MUST  BE  GREATER  THAN  OR  EQUAL  TO  ZERO 

SUBROUTINES  AND  FUNCTION  SUBPROGRAMS  REQUIRED 
NONE 

METHOD 

RECURRENCE  RELATION  AND  TRIG  FORMULAE  AS  DESCRIBED 
BY  H.A.  ANTQSIEWICZ  IN  HANDBOOK  OF  MATHEMATICAL 
FUNCTIONS »  PP.  438-439.  THIS  CODE  PATTERNED  AFTER 
THAT  FOR  BESY. 


SUBROUTINE  SBESY (X »N.BY » 1ER) 

CHECK  FOR  ERRORS  IN  N  AND  X 

IF<N)180, 10.10 
10  IER=0 

IF<X)190, 190.20 

EVACUATE  YO  AND  Yl. 

20  Y0=-C0S<X)/X 

Y1=-C0S  <X)/(X*X)-SIN(X)/X 

CHECK  IF  ONLY  YO  OR  Yl  IS  DESIRED 

90  IF <N-1 ) 100 » 100. 130 

RETURN  EITHER  YO  OR  Yl  AS  REQUIRED 

100  IF(N)UO. 120.110 
110  BY-Y1 

GO  TO  170 
120  BY=Y0 

GO  TO  170 

PERFORM  RECURRENCE  OPERATIONS  TO  FIND  YN(X> 
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0013 

130 

YA*YO 

0014 

YB-Y1 

0013 

K«1 

0016 

140 

T*FL0AT<7«X+1)/X 

0017 

YOTBYB-YA 

0018 

IF ( ABS( YC>-1 «0E36) 145 >145 >141 

0019 

141 

IER-3 

0020 

RETURN 

0021 

145 

K-K+l 

0022 

IF (K-N ) 150  > 160  r 150 

0023 

150 

YA-YB 

0024 

YB-YC 

002S 

60  TO  140 

0026 

160 

BY-YC 

0027 

170 

RETURN 

0028 

180 

IER-1 

0029 

RETURN 

0030 

190 

IER-2 

0031 

RETURN 

0032 

END 

* 
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Complex  Spherical  Bessel  Functions 


CSBSJ 

This  procedure  evaluates  the  spherical  Bessel  of  the  first  kind, jn(Z) „ 
for  complex  arguments  by  means  of  a  modified  continued  fraction  technique. 

This  technique  is  far  more  effective  than  upward  or  downward  recursion  methods 
The  spherical  Bessel  function  of  order  n  can  be  expressed  in  terms  of 
smaller  orders  as 


j.u) 

,(Z)  =  " 


*n-l 


(Z) 


3“TTZ)  3T^r 


n-2 


J-|(Z) 

37ZT 


*  j0(Z). 


(1) 


and  j Q ( Z )  =  sin  (Z)/Z. 

Each  of  the  ratios  in  (1)  is  found  by  a  modified  continued  fraction 
method.  Now 


jj<z)  _  ww 

by  definition.  From  Abramowitz  and  Stegiin2, 

^i(Zj  ,111 

3771)  a0"  al"  a2~ 

where  ak  =  2(i  +  k)/Z. 

The  inverse  of  this  ratio  can  be  rewritten  as 

j’i_i  (z)  n0*n1*n2*  .  .  . 

jyz) 

where  Nk  =  Ak  -  1/N|c_1  ,  NQ  =  AQ 


(2) 


(3) 


(4) 


»j 


k-1*  uo 


In  the  algorithm,  the  computation  of  (4)  is  continued  until  VDk-l 
-1  <  e,  where  e  is  the  allowable  error. 

This  program  was  written  by  David  Standley. 
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0001 


0002 

0003 

0004 

0003 

0004 

0007 

0008 

0009 

0010 

0011 

0013 

0014 

0013 

0014 

0017 

0018 

0019 

0021 

0023 

0024 

0025 

0024 

0027 

0029 

0030 

0031 

0032 

0034 


0035 

0034 

0037 


SUBROUTINE  CSBSJ(Xr  JJ.  Mr  IER > 

C*****************»****X****4«*t*«****>************************t*****t****** 

C  12-FEB-79  CHECK  FOR  XR-O.  BEFORE  ALOG  IS  TAKEN 
B-DEC-78  STOP  ITERATIONS  WHEN  (RATIO-1 XEPS 
24-0CT-78  OPTIMIZE  FOR  SPEED 
4-AUG-78 
D.  STANDLEY 


COMPUTE  COMPLEX  SPHERICAL  BESSEL  FUNCTIONS  BY  USE  OF  CONTINUED  FRACTIONS. 
SEE!  'A  METHOD  OF  COMPUTING  SPHERICAL  BESSEL  FUNCTIONS  OF  COMPLEX 
ARGUMENT  WITH  TABLES'  BY  W.  J.  LENTZ 

SUBROUTINE  CSBSJ(X.  JJ.  M.  IER) 

DESCRIPTION  OF  PARAMETERS 

X  -THE  ARGUMENT  OF  THE  J  SPHERICAL  BESSEL  FUNCTION  DESIRED 
X  IS  COMPLEX 
JJ  -THE  RESULTANT  J 

JJ  IS  THE  ARRAY  OF 
JJ  IS  COMPLEX 

M  -THE  ORDER  OF  THE  BESSEL  FUNCTION  DESIRED 
IER  -ERROR  CODE  WHERE: 

IER-0  NO  ERROR 

IER  .LT.  0  UNDERFLOW  OCCURRED  AT  ORDER  (-IER) 

IER*2  IMAGINARY  PART  OF  X  IS  >80.  OR  <-80. 

IER-3  ORDER  >  49 

THE  ACCURACY  IS  3  TO  4  SIGNIFICANT  FIGURES 


SPHERICAL  BESSEL  FUNCTION 
M  ORDERS 


C*****«***ts*y********s*«x**tt***********t****«***s***xt*****«t»««»«**tt**«*» 

COMPLEX  X.  At  SBESt  OUO<30) »  NUNN.  NUMNM 1 »  DENN.  DENNN1 r  JJ(SO) 

COMPLEX  ONE .  TWO  > TWODX . CSX . R  A  T 1 0 

A<N>  -  TWODX  *  CMPLX(V4FL0AT (N-l ) .  ZERO) 

Cl  ■  1. 

EPSR  -  l.E-5 
EPSI  -  l.E-3 
ZERO  -  0. 

ONE  -  (1..0.  ) 

TWO  -  (2..0.) 

IF<M  .GT.  49)G0  TO  314  I  ORDER  TOO  BIG 

XR  «  REAL ( X ) 

XI  «  AIMAG(X) 

CSX  -  CSIN(X) 

VV  <■  FLOAT (M) 

IXR  ■  0 
ICSR  -  O 

IF ( ABS ( XI )  . GT .  80. )G0  TO  513  ! OVERFLOW  CONDITION! 

IF (XR  .Ed.  ZERO) GO  TO  4 

ICSR  »  IFIX(ALOG10< ABS (REAL (CSX ) )  )  > 

IXR  -  IFIX ( ALOG 1 0 ( ABS ( XR ) ) > 

4  ICSI  *  0 
IXI  -  0 

IF (XI  .EQ.  ZERO  > GO  TO  3 

ICSI  -  IFIX(AL0G10(ABS( AI MAG (CSX ) ) ) ) 

IXI  >  IFIX (ALOGIO (ABS (XI>)> 

5  V  =  W  4  1.5 

IF(VV  .EQ.  ZERO ) GO  TO  501  (CALCULATE  DIRECTLY 

TWODX  -  TWO  /  X 

C _ 

C  CALCULATE  THE  RATIO  OF  J(V-1)/J(V> 

C 

DO  500  IV  ■  1.  IFIX(VV)  ! GET  NECESSARY  RATIOS 

V  *  V  -  1. 

L  *  IFIX(V  -  0.5) 
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0038  NUHNM1  «  A(l> 

0039  DENNM1  »A<2) 

0040  QUO(L)  -  NUHNN1 

0041  N  ■  1 

0042  10  N  «  N  +  1 

0043  NUNN  *  A (N)  -  ONE  /  NUHNM1 

0044  OENN  -  A(N+1>  -  ONE  /  DENNH1 

0043  RATIO  -  NUNN  /  DENNM1 

0044  IF ( ABS(REAL (RATIO ) -Cl )  .GT.  EPSR)GO  TO  11  'CHECK  FOR  CONVERGENCE 

0048  IF <ABS(AIMAG( RATIO))  .LT.  EPSDGO  TO  500  'STOP  THE  FRACTION 

0030  11  QUO(L)  *  QUO < L >  *  RATIO  ! CONTINUE  THE  FRACTION 

0031  NUHNM1  «  NUNN 

0032  0ENNN1  -  DENN 

0033  GO  TO  10 

0034  300  CONTINUE 

C - 

C  CALCULATE  EACH  ORDER 
C 


0055 

501 

DO  400  I»l>  N+l 

0054 

SBES  -  ONE 

0037 

N  -  I  -  1 

0058 

IF ( N  .EQ.  0>G0  TO 

310 

0040 

DO  302  J>  It  N 

0061 

502 

SBES  »  SBES  *  ONE 

/  QUO ( J ) 

C _ 

C  FROH  HERE  TO  310  CHECK  FOR  UNDERFLOW 
C 

0042  IF (REAL ( SBES)  .EG.  ZERO)GO  TO  312 

0044  ISR  -  IFIX< ALOGIO <ABS<REAL< SBES) ) > ) 

0043  ISI  -  0 

0044  IF(XI  .EQ.  ZERO > GO  TO  304 

0048  303  ISI  -  IFIX ( ALOGIO ( ABS ( AINAG ( SBES ) ) > ) 

0049  304  IF((1SR+ICSR).GT.  36. AND. < ISR+ICSR-IXR) . GT .-36 )G0  TO  303 

0071  GO  TO  512  IREAL  UNDERFLOW 

0072  305  IF((ISI+ICSI).GT .-34. AND «  < ISI+ICSI-IXI ) .GT .-34) GO  TO  310 

0074  GO  TO  312  ! INAG.  UNDERFLOW 

C _ 

C  HERE  CONES  THE  ANSWER 
C 


0075 

510 

SBES  -  SBES  *  CSX  /  X 

0074 

400 

JJ(I)  -  SBES 

0077 

IER  -  0 

! NORNAL  RETURN 

0078 

RETURN 

0079 

312 

IER  »  -  N 

! UNDERFLOW  AT  ORDER  N 

0080 

RETURN 

0081 

313 

IER  -  2 

! OVERFLOW 

0082 

RETURN 

0083 

514 

IER  »  3 

! ORDER  TOO  BIG 

0084 

END 

* 
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Physical  Notes 


The  reflectivity  of  a  fluid  sphere  generally  increases  as  its  density 
and  sound  speed  contrasts  increase.  The  reflectivity  increases  most  dramatically 
when  the  values  for  the  sphere  are  less  than  those  for  the  medium. 

For  the  elastic  sphere,  a  shear  speed  which  is  very  small  with  respect 
to  the  compressional  speed  of  the  sphere  and  the  medium  will  lead  to  a 
reflectivity  which  is  indistinguishable  from  the  corresponding  fluid  case. 
Apparently  the  impedance  contrast  in  this  case  is  such  that  little  or  no 
energy  is  converted  to  shear. 

The  elastic  sphere  program  produces  strange  results  for  targets  with 
high  shear  speed.  Many  of  these  cases  can  be  ruled  out  physically  since 
the  shear  speed  of  a  substance  must  be  less  than  0.7  x  its  compressional 
speed.  For  cases  with  moderate  shear  speeds,  the  results  are  similar  to  those 
for  a  fluid  sphere  but  spikey.  It  is  probably  necessary  to  include  the 
effects  of  attenuation  in  order  to  get  valid  results. 

For  cases  with  very  small  attenuation,  the  viscoelastic  sphere  program 
produces  results  which  are  indistinguishable  from  the  corresponding  elastic 
cases.  A  moderate  amount  of  shear  attenuation  will  smooth  out  the  spikiness 
caused  by  a  moderate  shear  speed.  Moderate  values  of  compressional  attenuation 
will  raise  the  level  of  the  curve;  large  values  will  also  decrease  the  depth 
of  the  dips. 
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Computational  Notes 


These  programs  and  subroutines  have  explicit  parameters  that  are  used 
to  test  for  convergence.  In  order  to  reduce  execution  times,  these  parameters 
have  been  set  relatively  high.  The  values  were  selected  on  the  basis  of  the 
output  plots,  so  the  relative  accuracy  may  be  only  about  one  percent. 

The  output  plots  are  generated  by  means  of  an  in-house  developed 
plotting  package.  The  output  sections  of  the  programs  can  be  easily  mod¬ 
ified  to  work  with  other  locally  available  plotting  software. 
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